#!/usr/bin/env julia
"""
Diagnostic File Generation Demo

This demo shows how to generate GSI-compatible diagnostic files using
the built-in Diagnostics module from GSICoreAnalysis.jl.
"""

using Printf
using Dates
using Random

using GSICoreAnalysis

const Diagnostics = GSICoreAnalysis.Diagnostics

"""
    generate_synthetic_observations(n_obs::Int, grid_size::Tuple)

Generate synthetic observation data for testing.
"""
function generate_synthetic_observations(n_obs::Int, grid_size::Tuple)
    println("Generating $n_obs synthetic observations...")

    Random.seed!(42)  # Reproducible

    observations = Diagnostics.ObservationRecord[]

    for i in 1:n_obs
        obs = Diagnostics.ObservationRecord()

        # Random location
        obs.lat = rand() * 180.0 - 90.0  # -90 to 90
        obs.lon = rand() * 360.0         # 0 to 360
        obs.pressure = rand() * 1000.0 + 10.0  # 10 to 1010 hPa

        # Synthetic atmospheric values
        obs.obs_value = 1013.0 + randn() * 5.0  # Surface pressure in hPa
        obs.obs_error = 1.0 + rand() * 0.5

        # Background and analysis values
        obs.background_value = obs.obs_value + randn() * 2.0
        obs.analysis_value = obs.obs_value + randn() * 1.0

        # Compute innovations
        obs.ob_minus_background = obs.obs_value - obs.background_value
        obs.ob_minus_analysis = obs.obs_value - obs.analysis_value

        # QC and flags
        obs.qc_flag = rand() < 0.95 ? Int32(0) : Int32(1)  # 5% rejected
        obs.use_flag = obs.qc_flag == 0 ? Int32(1) : Int32(0)

        # Observation type
        obs.obs_type = Diagnostics.get_obs_type_code("ps")
        obs.obs_subtype = Diagnostics.get_obs_subtype_code("adpsfc")

        # Station ID
        obs.station_id = @sprintf("STN%05d", i)

        # Time offset
        obs.time_offset = randn() * 0.5  # ±0.5 hours

        # Inverse error and weight
        obs.inverse_obs_error = 1.0 / obs.obs_error
        obs.variational_qc_weight = obs.use_flag > 0 ? 1.0 : 0.0

        push!(observations, obs)
    end

    println("  Generated $(length(observations)) observations")
    println("  Used: $(count(obs -> obs.use_flag > 0, observations))")
    println("  Rejected: $(count(obs -> obs.qc_flag > 0, observations))")

    return observations
end

"""
    generate_diagnostic_files_demo()

Main demo function to generate diagnostic files.
"""
function generate_diagnostic_files_demo()
    println("\n" * "="^80)
    println("GSI DIAGNOSTIC FILE GENERATION DEMO")
    println("="^80)
    println("Julia GSI Implementation")
    println("Time: $(Dates.format(now(), "yyyy-mm-dd HH:MM:SS"))")
    println("="^80 * "\n")

    # Configuration
    output_dir = joinpath(@__DIR__, "..", "results", "diagnostics")
    mkpath(output_dir)

    grid_size = (95, 57, 32)  # Match Fortran GSI example
    n_observations = 9000      # ~9000 observations like Fortran
    analysis_time = DateTime(2018, 8, 12, 12, 0, 0)

    println("Configuration:")
    println("  Output directory: $output_dir")
    println("  Grid size: $grid_size")
    println("  Number of observations: $n_observations")
    println("  Analysis time: $analysis_time")
    println()

    # Create diagnostic output container
    diag = Diagnostics.DiagnosticOutput(output_dir, analysis_time, grid_size)

    # Generate synthetic observations
    observations = generate_synthetic_observations(n_observations, grid_size)

    # Add observations to diagnostic container
    diag.observations["ps"] = observations

    # Compute innovation statistics
    stats = Diagnostics.compute_innovation_statistics(observations)
    diag.innovation_stats["ps"] = stats

    println("\nInnovation Statistics:")
    println("  Mean O-B: $(@sprintf("%.4f", stats["mean_innovation"]))")
    println("  RMS O-B:  $(@sprintf("%.4f", stats["rms_innovation"]))")
    println("  Bias:     $(@sprintf("%.4f", stats["bias"]))")
    println("  Std Dev:  $(@sprintf("%.4f", stats["std_innovation"]))")

    # Compute QC statistics
    diag.qc_stats["ps"] = Dict(
        "total" => length(observations),
        "used" => stats["n_used"],
        "rejected" => stats["n_rejected"]
    )

    # Add convergence info (synthetic)
    diag.iterations = 10
    diag.converged = true
    diag.cost_history = [100.0, 75.0, 55.0, 42.0, 35.0, 31.0, 29.0, 28.0, 27.5, 27.2, 27.1]
    diag.gradient_norms = [10.0, 7.5, 5.0, 3.5, 2.0, 1.5, 1.0, 0.7, 0.5, 0.3, 0.2]

    # Write all diagnostic files
    println("\n" * "="^80)
    println("WRITING DIAGNOSTIC FILES")
    println("="^80)

    Diagnostics.write_all_diagnostics(diag)

    # List generated files
    println("\n" * "="^80)
    println("GENERATED FILES")
    println("="^80)

    files = readdir(output_dir)
    for file in sort(files)
        filepath = joinpath(output_dir, file)
        size = filesize(filepath)
        @printf("  %-40s %10d bytes\n", file, size)
    end

    println("\n" * "="^80)
    println("✓ DEMO COMPLETE")
    println("="^80)
    println("\nDiagnostic files written to: $output_dir")
    println("\nTo validate against Fortran GSI, run:")
    println("  julia validation/validate_diagnostic_files.jl")
    println("="^80 * "\n")
end

# Run demo
if abspath(PROGRAM_FILE) == @__FILE__
    generate_diagnostic_files_demo()
end
